% Function : Forward Search Algorithm with Correlation
% Parameter:
%   Y : target data that we are trying to predict
%   X : training set, each row is one training example, each column is one feature
%   u : upper limit of selected features, search depth
%   gc: values of gray correlation
% Return:
%   I: indexes of selected features
%   R : average residuals of selected features
%   RR: R^2 statistic of selected features
function [I, R, RR] = CForward(Y, X, u, gc)
    if nargin ~= 4
        error('Usage: Forward(Y, X, u, gc)');
    end

    [m, ~] = size(X);   % m: number of samples, n: number of features

    [~, L] = KMeans(gc, 3);

    % initialize collection
    F  = ones(m, u + 1);    % collection of features
    I  = zeros(1, u);       % index of collected features
    R  = zeros(1, u);       % average residuals
    RR = zeros(1, u);       % R^2 statistic

    i = 0;
    k = 1;
    while k <= u
        i = mod(i, 3) + 1;  % i'th cluster

        % indexes of features in i'th cluster
        ft  = find(L == i);
        nft = numel(ft);    % number of features

        r  = zeros(1, nft);
        rr = zeros(1, nft);

        % traverse each feature in i'th cluster
        for j = 1 : nft
            % feature ft(j) is already in collection?
            idx = ft(j);
            if any(idx == I)
                continue;
            end

            F(:, k + 1) = X(:, idx);
            [~, ~, t, ~, stats] = regress(Y, F(:, 1 : k + 1), 0.05);
            r(j)  = mean(t .^ 2);
            rr(j) = stats(1);
        end

        if all(rr == 0)
            continue;
        end

        % select the feature that has the max R^2 statistic
        maxi = find(rr == max(rr));
        maxi = maxi(1);

        % udpate collection
        F(:, k + 1) = X(:, ft(maxi));
        I(k)  = ft(maxi);
        R(k)  = r(maxi);
        RR(k) = rr(maxi);

        k = k + 1;
    end
end
